#this is a script for customer spider/radar plot adapted from https://github.com/region-spotteR/PrepPlot
#library(plotly)
library(dplyr)
library(ggplot2)
library(raster)
library(data.table)


#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"


#===================Input data==================
load(paste0(wd1,"/data/RData/Fig.ED1-2.Radar-regional-effect.RData"))

eff.tot.all[Region!="Global"]
eff.tot.all[Region=="Global"] #global total effects

#===option 1: only show climate+CO2.eff for ER (RCP45(SSP5LU)- RCP85)
reg.tot <- copy(eff.tot.all[Region!="Global"])
#include co2 effect for ER scenario
reg.tot <- reg.tot[, Y.mod.tot:=(yield-yield0)/yield0*100] #relative change of global average yield (yield=RCP45_85LU,SAI,MSB,CCT, yield0=RCP85)
df.reg <- eff.crop.all[Region!="Global", .(Scenario,Crop,Region,T.eff, R.eff, M.eff, co2.eff)]

# #===option 2: show both CO2.eff and luc.eff for ER 
# reg.tot <- copy(eff.tot.all[Region!="Global"])
# reg.tot <- reg.tot[Scenario=="RCP45 - RCP85", Y.mod.tot:=Y.mod.dif3] #include co2.eff and luc.eff
# reg.tot <- reg.tot[Scenario!="RCP45 - RCP85", Y.mod.tot:=Y.mod.dif]
# df.reg <- eff.crop.all[Region!="Global", .(Scenario,Crop,Region,T.eff, R.eff, M.eff, co2.eff, luc.eff, Y.mod.dif)]

df.area <- eff.crop.all[Scenario=="RCP45 - RCP85", c("Crop", "area", "Region")]
setorder(df.area, Region, -area) #descending order by area
df.area <- df.area[,order:=order(area,decreasing=T), by=.(Region)]
print(df.area)

#pick two typical crops per region to get 12 region-crop subsets for radar plot:
df.reg <- df.reg[Crop=="Corn", Crop:="Maize"]
df.reg <- df.reg[Crop=="Soybean", Crop:="Soy"]
df.reg <- df.reg[Crop=="Sugarcane", Crop:="Cane"]
#pick 3 typical crops per region 
df.radar <- rbind(df.reg[Region=="NAM" & Crop%in%c("Wheat", "Maize", "Soy")],
                  df.reg[Region=="CAM" & Crop%in%c("Maize", "Cane", "Rice")],
                  df.reg[Region=="SAM" & Crop%in%c("Maize", "Soy", "Cane")], #Rice slightly higher than sugarcane
                  df.reg[Region=="SAF" & Crop%in%c("Maize", "Rice", "Cotton")], #Wheat higher than cotton
                  df.reg[Region=="EUA" & Crop%in%c("Wheat", "Maize", "Cotton")],
                  df.reg[Region=="SEA" & Crop%in%c("Rice", "Wheat", "Maize")])
scenario <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85")
#name <- c("Emissions reduction", "SAI - RCP85", "MSB - RCP85", "CCT - RCP85")
name <- c(expression(''*bold(a)*'     SAI - RCP85'), expression(''*bold(b)*'      MSB - RCP85'), 
          expression(''*bold(c)*'     CCT - RCP85'), expression(''*bold(d)*'      Emissions reduction'))
df.radar <- df.radar[reg.tot, Reg.tot:=i.Y.mod.tot, by=.EACHI, on=.(Scenario,Region)] #add regional tot effect
df.radar <- df.radar %>% 
  mutate(Region = factor(Region, levels=c("NAM", "CAM", "SAM", "SAF", "EUA", "SEA"))) %>%
  mutate(Scenario = factor(Scenario, levels=scenario, labels=name )) %>% as.data.table()
#make variable into a factor with the level in self-defined order and then setorder()
setorder(df.radar, Scenario, Region, Crop)
df.radar$Reg.tot[-seq(2,71,3)] <- NA #show only mid point per region by geom_path


#combine region-crop as group
df.radar$RegCrop <- paste(as.character(df.radar$Region),c(df.radar$Crop), sep="-") 
#===option 1: only include co2.eff for ER
df.radar2 <- df.radar[,c("Scenario", "RegCrop","T.eff","R.eff","M.eff","co2.eff", "Reg.tot")]
#===option 2: show both CO2.eff and luc.eff for ER
#df.radar2 <- df.radar[,c("Scenario", "RegCrop","T.eff","R.eff","M.eff","co2.eff","luc.eff", "Reg.tot")]

#use 12 RegCrop as axises instead of climate variables
df.radar.t <- df.radar2[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"RegCrop", by=.(Scenario)]
names(df.radar.t)[-1] <- c("Variable", unique(df.radar2$RegCrop))


#====plot all 4 scenarios in one panel
df_radar.all <- df.radar.t%>%
  as_tibble() %>% 
  group_by(Scenario,Variable) 

### 0.4 Set the plot parameters - mostly similar to ggradar
plot.data <- df_radar.all 
colnames(plot.data)
group_id <- c(1:2) 
group_n <- length(group_id)

(var.names <- colnames(plot.data)[-group_id])
(nvar <- length(var.names))
#axis.labels=c("'Temperature'", "'Radiation'", "'Moisture'", "'CO'['2']*''", "Land-use", "Total effect")
(axis.labels=colnames(plot.data)[-group_id])
(reg.labels <- substr(axis.labels, 1,3)) #take 1:3 characters
(crop.labels <- substring(axis.labels, 5)) #take all characters from 5 to end
#var.labels=c("Temperature effect", "Radiation effect", "Moisture effect", expression('CO'[2]~' effect'), "Land-use effect", "Regional total effect")
var.labels=c("Temperature effect", "Radiation effect", "Moisture effect", expression('CO'[2]~' effect'), "Regional total effect")



#####################################################################################################
############################ 0. Helper functions ####################################################
#####################################################################################################

#x y coordinates for our data, each angle for a variable
CalculateGroupPath4 <- function(df) {
  angles = seq(from=0, to=2*pi, by=(2*pi)/nvar) ##nvar+1 angles
  xx<-rbind(t(df[,-group_id])*sin(angles[-length(angles)]),t(df[,group_n+1])*sin(angles[1])) #nvar+1 lines for each group
  yy<-rbind(t(df[,-group_id])*cos(angles[-length(angles)]),t(df[,group_n+1])*cos(angles[1])) #each column for a group
  group <- df[,group_id] #array same ncol as xx,yy
  group_df <- as.data.frame(group[rep(seq_len(nrow(group)), each = length(angles)),]) 
  #if (group_n==1) graphData<-data.frame(group=rep(df[[group_id]],each=ncol(df)),x=c(xx),y=c(yy)) else
  graphData<- cbind(group_df, data.frame(x=c(xx),y=c(yy)))
  return(graphData)
}

#each axis corresponding to a variable (partial effects)
CalculateAxisPath2 <- function(var.names,min,max) {
  n<-length(var.names)
  #Cacluate required number of angles (in radians)
  angles <- seq(from=0, to=2*pi, by=(2*pi)/n)
  #calculate vectors of min and max x+y coords
  min.x <- min*sin(angles)
  min.y <- min*cos(angles)
  max.x <- max*sin(angles)
  max.y <- max*cos(angles)
  tmp<-lapply(1:n,function(i) matrix(c(i,i,min.x[i],max.x[i],min.y[i],max.y[i]),2,3))
  res<-as.data.frame(do.call(rbind,tmp))
  colnames(res) <- c("axis.no","x","y")
  return(res)
}

funcCircleCoords <- function(center = centre.y, r = 1, npoints = nvar+1 ){ #npoints = ncol(plot.data)){
  #Adapted from Joran's response to http://stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2
  tt <- seq(0,2*pi,length.out = npoints)
  yy <- center + r * cos(tt)
  xx <- center + r * sin(tt)
  return(data.frame(x = xx, y = yy))
}


#####################################################################################################
############################ 1. set common parameters ###############################################
#####################################################################################################

grid.min=-50 # 0  
grid.mid=0 #0.5
grid.max=100 #1 
#centre.y=grid.min - ((1/9)*(grid.max-grid.min))
centre.y=grid.min - ((1/6)*(grid.max-grid.min))
plot.extent.x.sf=1.25 #1.2  
plot.extent.y.sf=1.25
axis.label.offset= 1.1 #0.9
grid.label.offset= .9 #1.1
axis.line.colour="grey"
grid.line.colour="darkgrey"
background.circle.transparency=0.4

#===Radius of the gridlines
#r<-seq(0,1,0.25) 
r0 <- seq(grid.min,grid.max,50) #corresponding to the grid min, mid, max
r0 <- c(r0, grid.max*1.3) #add a outer band for labeling six regions


#####################################################################################################
############################ 2. Prepare Plotting  ####################################################
#####################################################################################################

### 2.1. Create vector with all KPI/variable names; Set up the data for plotting 
#var.names <- colnames(plot.data)[-1]  # Short version of variable names 
plot.data.offset <- plot.data
plot.data.offset[,-group_id]<- plot.data[,-group_id]+abs(centre.y)

### 2.2. Calculate the x and y coordinates for our data
xy_lines <- CalculateGroupPath4(plot.data.offset)
#xy_lines can be a dataframe with >1 group variables for faceting plot
if (group_n>1) xy_lines %>% 
  mutate_at(vars(Scenario), funs(factor(., levels=unique(.)))) %>%
  mutate_at(vars(Variable), funs(factor(., levels=unique(.))))-> xy_lines # convert to factor
  
### 2.3. Create a list containing all grid-objects:               
## 2.3.1 Calculate the data frame for the axis-lines (each axis corresponds to a crop)
grid <- NULL
grid$axis_path  <- CalculateAxisPath2(var.names,grid.min+abs(centre.y),grid.max+abs(centre.y))
if (group_n>1) {  #add scenario as group for facet plot
  n.r <- nrow(grid$axis_path)
  grid$axis_path <- grid$axis_path[rep(seq_len(n.r), time=length(unique(plot.data$Scenario))), ]
  grid$axis_path$Scenario <- factor(rep(unique(plot.data$Scenario), each=n.r), levels=unique(plot.data$Scenario)) 
}

## 2.3.2 Calculate the coordinates for the axis labels
grid$axis_label <- funcCircleCoords(0,(grid.max+abs(centre.y))*axis.label.offset,(nvar+1))[-(nvar+1),]
grid$axis_label$text=axis.labels  #use special characters later

## 2.3.3b For circular radar: Calculate the grid-lines (3 or 5 circles)
grid$lines_circle<- lapply(1:length(r0),function(i) funcCircleCoords(0,r0[i]+abs(centre.y),(nvar+1)*(2^7)))
names(grid$lines_circle)<- paste(round(r0,0),'%',sep='')
#names(grid$lines_circle)<-paste(seq(grid.min,grid.max,25),'%',sep='')

#2.3.4 Calculate the coordinates for grid circle labels
grid_label <- lapply(1:length(r0),function(i) funcCircleCoords(0,(r0[i]+abs(centre.y))*grid.label.offset,36)[2,]) #9 devides a circle to 9 sections each 45 degree, and take the 2nd row for points at 45 degree instead of labeling at 90 degree (x=0)
grid$grid_label <- do.call(rbind.data.frame, grid_label) #turn list to data.frame
grid$grid_label$text= paste(round(r0, 0),'%',sep='')
#if (group_n>1) grid$grid_label$Scenario <- factor(rep(unique(plot.data$Scenario), each=nrow(r2)),levels=unique(plot.data$Scenario))  
#not needed

## 2.3.5b For a circular grid: Bind all the grid-lines in one df and add the real values
grid$all_lines_c<-do.call(rbind, grid$lines_circle)
n<-nrow(grid$all_lines_c)/length(grid$lines_circle)
grid$all_lines_c$q<-rep(names(grid$lines_circle),each=n) # The % range of each grid
#if (group_n>1) grid$all_lines_c$Scenario <- factor(rep(unique(plot.data$Scenario), each=n*nrow(r2)),levels=unique(plot.data$Scenario))


#####################################################################################################
############################ 3. Plotting radar (common frame) ####################################################
#####################################################################################################

### Combine with all the code in the first chunk of shiny_example.Rmd to make it work
### 3.1. Create an empty theme for ggplot
theme_clear <- theme_bw(base_size=20) + 
  theme(axis.text.y=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.border=element_blank(),
        legend.key=element_rect(linetype="blank"))

### 3.2. Set the extent of the plot and the colors of the grid
plot.extent.x=(grid.max+abs(centre.y))*plot.extent.x.sf
plot.extent.y=(grid.max+abs(centre.y))*plot.extent.y.sf
#bg_colors<-c(RColorBrewer::brewer.pal(9,'Purples')[2],'#DFDFED') # Set the grid colors
bg_colors<- alpha(c("darkgreen"), background.circle.transparency)# c('gray95')

### 3.3. Set up a basic empty gg-object
basething <- ggplot() + xlab(NULL) + ylab(NULL) + coord_fixed() +
  scale_x_continuous(limits=c(-plot.extent.x,plot.extent.x)) + 
  scale_y_continuous(limits=c(-plot.extent.y,plot.extent.y)) + 
  theme_clear 

### 3.4 For the circular radar: Add the gridlines
n<-nrow(grid$all_lines_c)/length(grid$lines_circle)
c<-length(unique(plot.data$Scenario))
outer_circle <- grid$lines_circle[[length(r0)]]
nrow(outer_circle)/6/3 #135 points per crop segment; before divide the circle, shift half crop segment to the begining of point dataframe
outer_circle <- rbind(outer_circle[(nrow(outer_circle)-67):nrow(outer_circle),], outer_circle[1:(nrow(outer_circle)-68),])
div <- rep(seq_len(6), each = round(nrow(outer_circle)/6),length.out = nrow(outer_circle)) #factor divide the circle to six chuncks
outer_circle6 <- split(outer_circle, div)
#lapply(outer_circle6, dim)
outer_circle6$`1` <- outer_circle6$`1`[-c(406:407),] #remove the extra two points in list 1 
grid$outer_circle6 <- do.call(rbind.data.frame, outer_circle6) #turn list to data.frame
grid$outer_circle6$region= rep(c("NAM", "CAM", "SAM", "SAF", "EUA", "SEA"), as.integer(lapply(outer_circle6, nrow)))


library(RColorBrewer)
#colors <- brewer.pal(n =6, "Paired")
#reg.color <- rep(rep(c(colors[1],colors[c(4,6)],colors[2],colors[c(3,5)]), each=405),time=c) 
colors <- brewer.pal(n =8, "Pastel1")
reg.color <- rep(rep(colors[2:7], each=405),time=c)

base_grid_c <- basething + 
  geom_path(data=grid$outer_circle6,aes(x=x,y=y,group=region), color=reg.color, alpha=1,
            size=2,linetype=1, show.legend = F) + #draw color band for six regions
  #geom_polygon(data=grid$lines_circle[[length(r0)]],aes(x,y),fill="red3",linetype=1,colour="NA") + #draw outer color band for tot eff
  geom_polygon(data=grid$lines_circle$`100%`,aes(x,y),fill="white") + #make the inside hollow for the color band of tot.eff
  geom_polygon(#data=grid$all_lines_c[nrow(grid$all_lines_c):1,],aes(x,y,group=rev(q)), #from outer to inner
               data=grid$all_lines_c, aes(x,y,group=q), #draw polygon from inside out
               size=rep(rep(c(1,1.,1,1,0), each=n),c), colour="darkgray", fill=NA, #important not to fill at this step otherwise some circles will be blocked by others
               #colour=rep(rep(c("darkgray","darkgray","black","darkgray","darkgray", NA), each=n),c), 
               linetype=rep(rep(c(2,2,2,2,0), each=n),c))+ 
               #fill=rep(c(rep(bg_colors,each=n), rep('white',3*n)),c)) + #2 colors, outer band-red for tot net effect, inner circel-white for partial effects  
               #fill=rep(c(rep("white",length(r0)*n), rep("NA",each=n)),c)) + #colored band for tot net effect (last item)
  geom_path(data=grid$axis_path,aes(x=x,y=y,group=axis.no), colour=axis.line.colour,alpha=0.6,size=1,linetype="solid") +
  facet_grid(~ Scenario) 

### 3.5 Add text and label to gg-object
#for 3x6 crop-region combination
angles <- seq(90, -250,length.out = 18)
angles[11:18] <- angles[11:18]+180  #flip text angles for the left
vjust <- rep(0.5, 18) #rep(-0.2, 18)  
vjust[6:14] <- rep(0.3,9) #rep(1.,9)
angles2 <- seq(0,-340,length.out = 18)
angles2[6:14] <- angles2[6:14]+180 
vjust2 <- rep(-1.4, 18) # rep(-2.5, 18)  
vjust2[6:14] <- rep(2.4,9) # rep(3.5,9)

reg.labels[-seq(2,17,3)] <- NA #only label mid-point for each region

base_text <- base_grid_c+ 
  geom_text(data=grid$axis_label,aes(x,y,label=rep(crop.labels,c)), parse = TRUE, #display expressions by parse = TRUE
            vjust=rep(vjust,c), angle=rep(angles2,c), size=6, fontface=1)+   #hjust=rep(c(.5,.5,.5,.5,.5,.5),c) #crop label
  geom_text(data=grid$axis_label,aes(x,y,label=rep(reg.labels,c)), parse = TRUE, #show region label
            vjust=rep(vjust2,c), angle=rep(angles2,c), size=7, fontface=1)+   #hjust=rep(c(.5,.5,.5,.5,.5,.5),c), 
  geom_text(data=grid$grid_label,aes(x,y,label=text),size=rep(c(rep(6,4),0),c), fontface=rep(c(1,2,1,1,0),c),
            vjust=rep(c(1,.5,.5,0,0),c),hjust=rep(c(.5,.5,.5,.5,0),c)) #for 3x6 crop-region



#================================================================
#=============all 4 plots together for SGs and ER================
#================================================================
varcolor2 <- c("#DC0000FF", "#FF7F00", "#56B4E9", "#00A087FF","#000000")#
barplot(height=rep(1,length(varcolor2)),col=varcolor2) #display colors


### 3.6 Add data lines to gg-object 
xy_lines <- as.data.table(xy_lines)

if (c>1) {
  xy_lines[Variable=="Reg.tot"][(1:4)*(nvar+1),] <- xy_lines[Variable=="Reg.tot"][2+(0:3)*(nvar+1),] #need a extra point to close the path 
} 
reg.labels[-seq(2,17,3)] <- NA #only label mid-point for each region


#(!use na.omit(xy_lines) so that geom_path can connect across missing points)
gg<-base_text+geom_path(data=na.omit(xy_lines), aes(x=x,y=y,group=Variable,color=Variable),
                        alpha=0.6, size=1) + #size=rep(c(1,1,1,1,1,1.5),each=19)) +
  geom_point(data=xy_lines, na.rm = T, aes(x=x,y=y,group=Variable,color=Variable, shape=Variable), 
             shape=rep(rep(c(19,19,19,19,8),each=19),4), alpha=0.6, size=rep(rep(c(2,2,2,2,3)*1.2,each=19),4)) +
  scale_colour_manual(name="", labels=var.labels, values = c(varcolor2)) +
  scale_shape_manual(values = c(19,19,19,19,8))+
  #facet_grid(~ Scenario) + 
  facet_wrap(~Scenario, nrow=2, labeller= label_parsed)+
  #theme(panel.spacing = unit(5, "lines")) +  #adjust facet spacing (not helpful, adjust output jpg size)
  theme(strip.background = element_blank(), strip.text = element_text(size = 20))+
  guides(colour=guide_legend(title=NULL, nrow=2, byrow=T,  keywidth = 2, keyheight=1, label.position = "right", label.hjust = 0,
                             override.aes = list(linetype=1, shape=c(19,19,19,19,8), alpha=0.6 )))+
  #shape =guide_legend(override.aes = list(size = 10, shape=c(19,19,19,19,4), alpha=0.6))) +
  theme(legend.text = element_text(size=20), legend.position = "bottom") #,legend.key.size = unit(10, "points")) 
gg
##major crops
#ggsave(paste0(wd1,"/figures/extended/Fig.Ext1.svg"), device ="svg", units = "in", width = 14, height = 14, dpi = 600)
ggsave(paste0(wd1,"/figures/extended/ED_Fig1.jpg"), device ="jpg", units = "in", width = 14, height = 14, dpi = 350)







##################Extended figure 2 for other minor crops#######################


#the other minor crops not included in the primary plot
df.radar <- rbind(df.reg[Region=="NAM" & Crop%in%c("Rice", "Cotton", "Cane")], #sugarcane in NAM is highly biased
                  df.reg[Region=="CAM" & !Crop%in%c("Maize", "Cane", "Rice")],
                  df.reg[Region=="SAM" & !Crop%in%c("Maize", "Soy", "Cane")],
                  df.reg[Region=="SAF" & !Crop%in%c("Maize", "Rice", "Cotton")],
                  df.reg[Region=="EUA" & Crop%in%c("Cotton", "Soy", "Rice")], #only five crops in EUA, reuse cotton
                  df.reg[Region=="SEA" & !Crop%in%c("Rice", "Wheat", "Maize")])
scenario <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85")
# name <- c(expression('('*bold(a)*') SAI - RCP85'), expression('('*bold(b)*') MSB - RCP85'), 
#           expression('('*bold(c)*') CCT - RCP85'), expression('('*bold(d)*') Emissions reduction'))
df.radar <- df.radar[reg.tot, Reg.tot:=i.Y.mod.tot, by=.EACHI, on=.(Scenario,Region)] #add regional tot effect
df.radar <- df.radar %>% 
  mutate(Region = factor(Region, levels=c("NAM", "CAM", "SAM", "SAF", "EUA", "SEA"))) %>%
  mutate(Scenario = factor(Scenario, levels=scenario, labels=name )) %>% as.data.table()
#make variable into a factor with the level in self-defined order and then setorder()
setorder(df.radar, Scenario, Region, Crop)
df.radar$Reg.tot[-seq(2,71,3)] <- NA #show only mid point per region by geom_path


#combine region-crop as group
df.radar$RegCrop <- paste(as.character(df.radar$Region),c(df.radar$Crop), sep="-") 
#===option 1: only include co2.eff for ER
df.radar2 <- df.radar[,c("Scenario", "RegCrop","T.eff","R.eff","M.eff","co2.eff", "Reg.tot")]

#use 12 RegCrop as axises instead of climate variables
df.radar.t <- df.radar2[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"RegCrop", by=.(Scenario)]
names(df.radar.t)[-1] <- c("Variable", unique(df.radar2$RegCrop))


#====plot all 4 scenarios in one panel
df_radar.all <- df.radar.t%>%
  as_tibble() %>% 
  group_by(Scenario,Variable) 

### 0.4 Set the plot parameters - mostly similar to ggradar
plot.data <- df_radar.all 
colnames(plot.data)
group_id <- c(1:2) 
group_n <- length(group_id)

(var.names <- colnames(plot.data)[-group_id])
(nvar <- length(var.names))
(axis.labels=colnames(plot.data)[-group_id])
(reg.labels <- substr(axis.labels, 1,3)) #take 1:3 characters
(crop.labels <- substring(axis.labels, 5)) #take all characters from 5 to end
var.labels=c("Temperature effect", "Radiation effect", "Moisture effect", expression('CO'[2]~' effect'), "Regional total effect")



### 2.1. Create vector with all KPI/variable names; Set up the data for plotting 
#var.names <- colnames(plot.data)[-1]  # Short version of variable names 
plot.data.offset <- plot.data
plot.data.offset[,-group_id]<- plot.data[,-group_id]+abs(centre.y)

### 2.2. Calculate the x and y coordinates for our data
xy_lines <- CalculateGroupPath4(plot.data.offset)
#xy_lines can be a dataframe with >1 group variables for faceting plot
if (group_n>1) xy_lines %>% 
  mutate_at(vars(Scenario), funs(factor(., levels=unique(.)))) %>%
  mutate_at(vars(Variable), funs(factor(., levels=unique(.))))-> xy_lines # convert to factor

### 2.3. Create a list containing all grid-objects:               
## 2.3.1 Calculate the data frame for the axis-lines (each axis corresponds to a crop)
grid <- NULL
grid$axis_path  <- CalculateAxisPath2(var.names,grid.min+abs(centre.y),grid.max+abs(centre.y))
if (group_n>1) {  #add scenario as group for facet plot
  n.r <- nrow(grid$axis_path)
  grid$axis_path <- grid$axis_path[rep(seq_len(n.r), time=length(unique(plot.data$Scenario))), ]
  grid$axis_path$Scenario <- factor(rep(unique(plot.data$Scenario), each=n.r), levels=unique(plot.data$Scenario)) 
}

## 2.3.2 Calculate the coordinates for the axis labels
grid$axis_label <- funcCircleCoords(0,(grid.max+abs(centre.y))*axis.label.offset,(nvar+1))[-(nvar+1),]
grid$axis_label$text=axis.labels  #use special characters later

## 2.3.3b For circular radar: Calculate the grid-lines (3 or 5 circles)
grid$lines_circle<- lapply(1:length(r0),function(i) funcCircleCoords(0,r0[i]+abs(centre.y),(nvar+1)*(2^7)))
names(grid$lines_circle)<- paste(round(r0,0),'%',sep='')
#names(grid$lines_circle)<-paste(seq(grid.min,grid.max,25),'%',sep='')

#2.3.4 Calculate the coordinates for grid circle labels
grid_label <- lapply(1:length(r0),function(i) funcCircleCoords(0,(r0[i]+abs(centre.y))*grid.label.offset,36)[2,]) #9 devides a circle to 9 sections each 45 degree, and take the 2nd row for points at 45 degree instead of labeling at 90 degree (x=0)
grid$grid_label <- do.call(rbind.data.frame, grid_label) #turn list to data.frame
grid$grid_label$text= paste(round(r0, 0),'%',sep='')
#if (group_n>1) grid$grid_label$Scenario <- factor(rep(unique(plot.data$Scenario), each=nrow(r2)),levels=unique(plot.data$Scenario))  
#not needed

## 2.3.5b For a circular grid: Bind all the grid-lines in one df and add the real values
grid$all_lines_c<-do.call(rbind, grid$lines_circle)
n<-nrow(grid$all_lines_c)/length(grid$lines_circle)
grid$all_lines_c$q<-rep(names(grid$lines_circle),each=n) # The % range of each grid
#if (group_n>1) grid$all_lines_c$Scenario <- factor(rep(unique(plot.data$Scenario), each=n*nrow(r2)),levels=unique(plot.data$Scenario))


### Combine with all the code in the first chunk of shiny_example.Rmd to make it work
### 3.1. Create an empty theme for ggplot
theme_clear <- theme_bw(base_size=20) + 
  theme(axis.text.y=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.border=element_blank(),
        legend.key=element_rect(linetype="blank"))

### 3.2. Set the extent of the plot and the colors of the grid
plot.extent.x=(grid.max+abs(centre.y))*plot.extent.x.sf
plot.extent.y=(grid.max+abs(centre.y))*plot.extent.y.sf
#bg_colors<-c(RColorBrewer::brewer.pal(9,'Purples')[2],'#DFDFED') # Set the grid colors
bg_colors<- alpha(c("darkgreen"), background.circle.transparency)# c('gray95')

### 3.3. Set up a basic empty gg-object
basething <- ggplot() + xlab(NULL) + ylab(NULL) + coord_fixed() +
  scale_x_continuous(limits=c(-plot.extent.x,plot.extent.x)) + 
  scale_y_continuous(limits=c(-plot.extent.y,plot.extent.y)) + 
  theme_clear 

### 3.4 For the circular radar: Add the gridlines
n<-nrow(grid$all_lines_c)/length(grid$lines_circle)
c<-length(unique(plot.data$Scenario))
outer_circle <- grid$lines_circle[[length(r0)]]
nrow(outer_circle)/6/3 #135 points per crop segment; before divide the circle, shift half crop segment to the begining of point dataframe
outer_circle <- rbind(outer_circle[(nrow(outer_circle)-67):nrow(outer_circle),], outer_circle[1:(nrow(outer_circle)-68),])
div <- rep(seq_len(6), each = round(nrow(outer_circle)/6),length.out = nrow(outer_circle)) #factor divide the circle to six chuncks
outer_circle6 <- split(outer_circle, div)
outer_circle6$`1` <- outer_circle6$`1`[-c(406:407),] #remove the extra two points in list 1 
grid$outer_circle6 <- do.call(rbind.data.frame, outer_circle6) #turn list to data.frame
grid$outer_circle6$region= rep(c("NAM", "CAM", "SAM", "SAF", "EUA", "SEA"), as.integer(lapply(outer_circle6, nrow)))

base_grid_c <- basething + 
  geom_path(data=grid$outer_circle6,aes(x=x,y=y,group=region), color=reg.color, alpha=1,
            size=2,linetype=1, show.legend = F) + #draw color band for six regions
  geom_polygon(data=grid$lines_circle$`100%`,aes(x,y),fill="white") + #make the inside hollow for the color band of tot.eff
  geom_polygon(#data=grid$all_lines_c[nrow(grid$all_lines_c):1,],aes(x,y,group=rev(q)), #from outer to inner
    data=grid$all_lines_c, aes(x,y,group=q), #draw polygon from inside out
    size=rep(rep(c(1,1.,1,1,0), each=n),c), colour="darkgray", fill=NA, #important not to fill at this step otherwise some circles will be blocked by others
    linetype=rep(rep(c(2,2,2,2,0), each=n),c))+ 
  geom_path(data=grid$axis_path,aes(x=x,y=y,group=axis.no), colour=axis.line.colour,alpha=0.6,size=1,linetype="solid") +
  facet_grid(~ Scenario) 

### 3.5 Add text and label to gg-object
#for 3x6 crop-region combination
angles <- seq(90, -250,length.out = 18)
angles[11:18] <- angles[11:18]+180  #flip text angles for the left
vjust <- rep(0.5, 18) #rep(-0.2, 18)  
vjust[6:14] <- rep(0.3,9) #rep(1.,9)
angles2 <- seq(0,-340,length.out = 18)
angles2[6:14] <- angles2[6:14]+180 
vjust2 <- rep(-1.4, 18) # rep(-2.5, 18)  
vjust2[6:14] <- rep(2.4,9) # rep(3.5,9)

reg.labels[-seq(2,17,3)] <- NA #only label mid-point for each region

base_text <- base_grid_c+ 
  geom_text(data=grid$axis_label,aes(x,y,label=rep(crop.labels,c)), parse = TRUE, #display expressions by parse = TRUE
            vjust=rep(vjust,c), angle=rep(angles2,c), size=6, fontface=1)+   #hjust=rep(c(.5,.5,.5,.5,.5,.5),c) #crop label
  geom_text(data=grid$axis_label,aes(x,y,label=rep(reg.labels,c)), parse = TRUE, #show region label
            vjust=rep(vjust2,c), angle=rep(angles2,c), size=7, fontface=1)+   #hjust=rep(c(.5,.5,.5,.5,.5,.5),c), 
  geom_text(data=grid$grid_label,aes(x,y,label=text),size=rep(c(rep(6,4),0),c), fontface=rep(c(1,2,1,1,0),c),
            vjust=rep(c(1,.5,.5,0,0),c),hjust=rep(c(.5,.5,.5,.5,0),c)) #for 3x6 crop-region

varcolor2 <- c("#DC0000FF", "#FF7F00", "#56B4E9", "#00A087FF","#000000")#
barplot(height=rep(1,length(varcolor2)),col=varcolor2) #display colors


### 3.6 Add data lines to gg-object 
xy_lines <- as.data.table(xy_lines)

if (c>1) {
  xy_lines[Variable=="Reg.tot"][(1:4)*(nvar+1),] <- xy_lines[Variable=="Reg.tot"][2+(0:3)*(nvar+1),] #need a extra point to close the path 
} 
reg.labels[-seq(2,17,3)] <- NA #only label mid-point for each region


#(!use na.omit(xy_lines) so that geom_path can connect across missing points)
gg<-base_text+geom_path(data=na.omit(xy_lines), aes(x=x,y=y,group=Variable,color=Variable),
                        alpha=0.6, size=1) + #size=rep(c(1,1,1,1,1,1.5),each=19)) +
  geom_point(data=xy_lines, na.rm = T, aes(x=x,y=y,group=Variable,color=Variable, shape=Variable), 
             shape=rep(rep(c(19,19,19,19,8),each=19),4), alpha=0.6, size=rep(rep(c(2,2,2,2,3)*1.2,each=19),4)) +
  scale_colour_manual(name="", labels=var.labels, values = c(varcolor2)) +
  scale_shape_manual(values = c(19,19,19,19,8))+
  #facet_grid(~ Scenario) + 
  facet_wrap(~Scenario, nrow=2, labeller= label_parsed)+
  #theme(panel.spacing = unit(5, "lines")) +  #adjust facet spacing (not helpful, adjust output jpg size)
  theme(strip.background = element_blank(), strip.text = element_text(size = 20))+
  guides(colour=guide_legend(title=NULL, nrow=2, byrow=T,  keywidth = 2, keyheight=1, label.position = "right", label.hjust = 0,
                             override.aes = list(linetype=1, shape=c(19,19,19,19,8), alpha=0.6 )))+
  #shape =guide_legend(override.aes = list(size = 10, shape=c(19,19,19,19,4), alpha=0.6))) +
  theme(legend.text = element_text(size=20), legend.position = "bottom") #,legend.key.size = unit(10, "points")) 
gg

#ggsave(paste0(wd1,"/figures/extended/Fig.Ext2.svg"), device ="svg", units = "in", width = 14, height = 14, dpi = 600)
ggsave(paste0(wd1,"/figures/extended/ED_Fig2.jpg"), device ="jpg", units = "in", width = 14, height = 14, dpi = 350)










